#Authors: Natili, Pellegata, Visconti
#Date: 06 November 2025
#Data: Eco-Social UNIMI (Ferrera et al. 2022)
#Paper Title: What Kind of Energy Transition? Public Opinion Trade-Offs among Economic
#       Growth, Ecological Sustainability and Equity
#Journal: European Journal of Political research
#Object: Create Figure 2 to 12 in the main text

#------------------------------------------------------------------------------#
#LOAD/INSTALL REQUIRED PACKAGES
#------------------------------------------------------------------------------#

#Load the relevant packages or install them if not available
if (!require("pacman")) install.packages("pacman")
pacman::p_load(cregg, tidyverse, ggplot2, patchwork)

#If cregg does not install please use the the following code:
#remotes::install_github("leeper/cregg")

#------------------------------------------------------------------------------#
#SPECIFY PATH TO WORKING DIRECTORY WITH DATASET
#------------------------------------------------------------------------------#
#Set working directory Windows
#setwd("C:\Users\fvisconti\Documents\Replication")

#Set working directory Mac OS
setwd("/Users/fv/Library/CloudStorage/Dropbox/ActiveWorks/myREScUE/Papers/Eco-social/eco-social-shared/Paper/EJPR-submission/R2/Replication")
#------------------------------------------------------------------------------#

#Clear all objects in environment
rm(list = ls())

#------------------------------------------------------------------------------#
#LOAD DATASET
#------------------------------------------------------------------------------#
#Load data object in long format ready for analysis
load(file = "dataset.RData")

#------------------------------------------------------------------------------#
#CREATE FIGURES FOLDER
#------------------------------------------------------------------------------#

dir.create("figures", showWarnings = FALSE)

#------------------------------------------------------------------------------#
##### EXCLUDE INCONSISTENT RESPONDENTS #####
#(chosed one package but rated the other one as favoured)
#------------------------------------------------------------------------------#
long <- long %>%
  group_by(pair) %>%
  mutate(inconsistent = ifelse(
    any(profile == 'A' & choice == 2 & rating > max(rating[profile == 'B' & choice == 2])) |
      any(profile == 'A' & choice == 1 & rating < max(rating[profile == 'B' & choice == 1])),
    1, 0
  )) %>%
  ungroup()

table(long$inconsistent)

#Drop inconsistent values
dim(long[long$inconsistent == 0,])

long <- long[long$inconsistent == 0,]

dim(long)

#------------------------------------------------------------------------------#
##### BASE MODEL - DISCRETE CHOICE WITH WEIGHTS #####
#------------------------------------------------------------------------------#

fit1mm <- chosen_rec ~ Energy*Compensation + Financing
mms <- mm(long, fit1mm,
          weights = ~weight, id = ~caseid)

################################################################################
# FIGURE 2 - MARGINAL MEANS
################################################################################

#Define font as Bold for the three factors (from bottom up)
myFaces <- c(rep('plain', 4), "bold",
             rep('plain', 4), "bold",
             rep('plain', 3), "bold")

#Figure
p1 <- plot(mms, vline = 0.5, header_fmt = "%s", size = 3,
           group = "feature") + 
  aes(shape = feature) + # map group to `shape` aesthetic
  scale_shape_manual(values=c(15, 16, 17)) +
  ggplot2::theme(legend.position = "none",
                 panel.grid.major = ggplot2::element_blank(),
                 panel.grid.minor = ggplot2::element_blank(),
                 axis.text.y = element_text(face=myFaces, size = 11)) +
  ggplot2::geom_errorbar(orientation= "y",
                          aes(xmin = lower,
                              xmax = upper),
                          linewidth = 0.75, height = 0, na.rm = TRUE) +
  labs(title = "Effects of attributes on support for a policy package",
       #       caption = "Source: Unimi Eco-Social Survey."
  )

p1
ggsave(filename = 'figures/Figure_2.png',
       plot = p1, width = 8, height = 6)


#------------------------------------------------------------------------------#
##### INTERACTIONS ####
#------------------------------------------------------------------------------#

#------------------------------------------------------------------------------#
# Interaction between energy and compensation
#------------------------------------------------------------------------------#

#Compute marginal means
interfit1 <- cj(long, chosen_rec ~ Compensation,
                id = ~ caseid,
                weights = ~ weight,
                estimate = "mm",
                h0 = 0.5,
                by = ~Energy)

#Define plot
interplot1 <- plot(interfit1, vline = 0.5,
                   header_fmt = "%s") +
  ggplot2::theme(plot.title = element_text(size=16),
                 plot.subtitle = element_text(size=12),
                 legend.position = "none",
                 axis.text=element_text(size=12),
                 axis.title=element_text(size=14),
                 axis.text.y = element_text(face=myFaces),
                 text = element_text(size = 12)) +
  geom_point(size=4) +
  labs(title = "Support for a policy package",
       subtitle = "Interaction between Energy and Compensation",
       #caption = "Source: Unimi Eco-Social Survey."
  )  + 
  ggplot2::facet_wrap(~BY, ncol = 4L) +
  ggplot2::geom_errorbarh(aes(xmin = lower,
                              xmax = upper),
                          linewidth = 0.75, height = 0, na.rm = TRUE)

interplot1
#Save plot
ggsave(filename = 'figures/Figure_3.png',
       plot = interplot1, width = 12, height = 6)

#------------------------------------------------------------------------------#
# Interaction between energy and financing
#------------------------------------------------------------------------------#

#Compute marginal means
interfit2 <- cj(long, chosen_rec ~ Financing,
                id = ~ caseid,
                weights = ~ weight,
                estimate = "mm",
                h0 = 0.5,
                by = ~Energy)

#Define plot
interplot2 <- plot(interfit2, vline = 0.5,
                   header_fmt = "%s") +
  ggplot2::theme(plot.title = element_text(size=16),
                 plot.subtitle = element_text(size=12),
                 legend.position = "none",
                 axis.text=element_text(size=12),
                 axis.title=element_text(size=14),
                 axis.text.y = element_text(face=myFaces),
                 text = element_text(size = 12)) +
  geom_point(size=4) +
  labs(title = "Support for a policy package",
       subtitle = "Interaction between Energy and Financing",
       #caption = "Source: Unimi Eco-Social Survey."
  )  + 
  ggplot2::facet_wrap(~BY, ncol = 4L) +
  ggplot2::geom_errorbarh(aes(xmin = lower,
                              xmax = upper),
                          linewidth = 0.75, height = 0, na.rm = TRUE)
interplot2
#Save plot
ggsave(filename = 'figures/Figure_4.png',
       plot = interplot2, width = 12, height = 6)

#------------------------------------------------------------------------------#
#### Occupation ####
#------------------------------------------------------------------------------#

#Compute marginal means
occ3 <- cj(
  subset(long, !is.na(long$occ2)),
  chosen_rec ~ Energy*Compensation + Financing, id = ~caseid, by = ~ occ2,
  weights = ~ weight, estimate = "mm"
)

#Exclude from plot Homemakers, Students, Not working/Others
occ4 <- occ3[!(occ3$BY %in% c("Homemaker",
                              "Not in work/Others",
                              "Student")),]

#Define style of text in figure
myFaces <- c(rep('plain', 4), "bold",
             rep('plain', 4), "bold",
             rep('plain', 3), "bold")

#Define the plot
occ2plot <- plot(occ4, group = "occ2",
                 vline = 0.5, 
                 header_fmt = "%s",
                 size = 3) +
  aes(shape = feature, group = "feature") + # map group to `shape` aesthetic
  scale_shape_manual(values=c(15, 16, 17)) +
  labs(title = "Support for policy package by occupation") +
  guides(colour=guide_legend(title="Occupation")) +
  ggplot2::theme(plot.title = element_text(size=20),
                 plot.subtitle = element_text(size=18),
                 axis.text=element_text(size=14),
                 axis.title=element_text(size=16),
                 axis.text.y = element_text(face=myFaces),
                 text = element_text(size = 18), 
                 legend.position = "none") +
  facet_wrap(~BY, ncol = 3L) +
  ggplot2::geom_errorbarh(aes(xmin = lower,
                              xmax = upper),
                          linewidth = 0.75, height = 0, na.rm = TRUE)

occ2plot

#Save the plot
ggsave(filename = 'figures/Figure_5.png',
       plot = occ2plot, width = 20, height = 12)


#------------------------------------------------------------------------------#
# Interaction between energy and compensation by occupation levels
#------------------------------------------------------------------------------#

# Initialize empty data frame
results_df <- data.frame()

# Loop over levels of occ2
for (occ_level in unique(long$occ2)) {
  # Subset data
  data_subset <- subset(long, occ2 == occ_level)
  
  # Run model
  model_output <- cj(data_subset, chosen_rec ~ Compensation,
                     id = ~ caseid,
                     weights = ~ weight,
                     estimate = "mm",
                     h0 = 0.5,
                     by = ~Energy)
  
  # Add column for occ2 group
  model_output$occ2 <- occ_level
  
  # Bind to results
  results_df <- bind_rows(results_df, model_output)
}


results_df2 <- results_df[!(results_df$occ2 %in% c("Homemaker",
                                                   "Not in work/Others",
                                                   "Student")),]


# Get unique values of occ2
occ2_vals <- unique(results_df2$occ2)

# Indexes where y-axis labels should be shown: 1st, 3rd, 5th
show_indices <- c(1, 3, 5)

# Create plot list
plot_list <- lapply(seq_along(occ2_vals), function(i) {
  o <- occ2_vals[i]
  show_y <- i %in% show_indices  # Only show y labels for specified positions
  
  ggplot(subset(results_df2, occ2 == o), aes(x = estimate, y = level)) +
    geom_point(size = 3, color = "red") +
    geom_errorbarh(aes(xmin = lower, xmax = upper),
                   linewidth = 0.75, height = 0, na.rm = TRUE,
                   color = "red")  +
    #geom_text(aes(label = sprintf("%.2f", estimate)),
    #          vjust = -1, size = 3, color = "black") +  
    facet_wrap(~BY, nrow = 1) +
    scale_x_continuous(limits = c(0.25, 0.75), breaks = seq(0.25, 0.75, by = 0.25)) +
    geom_vline(xintercept = 0.5, linetype = "dashed", color = "gray50") +
    labs(
      title = o,
      x = "Marginal Means",
      y = NULL
    ) +
    theme_bw(base_size = 12) +
    theme(
      plot.title = element_text(size = 10,face = "bold"),
      strip.text = element_text(size = 10),
      axis.text.x = element_text(size = 10),
      axis.title.x = element_text(size = 10),
      panel.grid = element_blank(),
      axis.text.y = if (show_y) element_text(size = 10) else element_blank(),
      axis.ticks.y = if (show_y) element_line() else element_blank()
    )
})

# Combine plots vertically
final_plot <- wrap_plots(plotlist = plot_list, ncol = 2) + 
  plot_annotation(
    title = "Interaction between energy and compensation by occupation levels"
  )
print(final_plot)

#Save plot
ggsave("figures/Figure_6.png",
       final_plot,
       width = 16, height = 2*length(occ2_vals))

#------------------------------------------------------------------------------#
#### Threatened by CC policies ####
#------------------------------------------------------------------------------#

#Compute marginal means
threatccpol1 <- cj(
  long,
  chosen_rec ~ Energy*Compensation + Financing, id = ~caseid, by = ~threatccpol,
  weights = ~ weight, estimate = "mm"
)

#Define plot
threatccpol1plot <- plot(threatccpol1[c(1:11, 23:33),], 
                         group = "threatccpol", 
                         header_fmt = "%s",
                         vline = 0.5, size = 3) +
  labs(title = "Support for policy package by job limited by CC policies") +
  guides(colour=guide_legend(title="Job limited by CC policies")) +
  aes(shape = feature) + # map group to `shape` aesthetic
  scale_shape_manual(values=c(15, 16, 17)) +
  ggplot2::theme(plot.title = element_text(size=20),
                 plot.subtitle = element_text(size=16),
                 axis.text=element_text(size=16),
                 axis.title=element_text(size=16),
                 text = element_text(size = 18), 
                 legend.position = "bottom",
                 legend.title = element_text(size=14),
                 legend.text = element_text(size=14),
                 axis.text.y = element_text(face=myFaces)
  ) +
  scale_shape(guide = 'none')
  

threatccpol1plot 

#Save plot
ggsave(filename = 'figures/Figure_7.png',
       plot = threatccpol1plot, width = 12, height = 10)


#------------------------------------------------------------------------------#
# Iteraction between energy, compensation and threatened by CC policiesFigure 8#
#------------------------------------------------------------------------------#

# Initialize empty data frame
results_df <- data.frame()

# Loop over levels of threatccpol
for (threatccpol_level in unique(long$threatccpol)) {
  # Subset data
  data_subset <- subset(long, threatccpol == threatccpol_level)
  
  
  # Run model safely
  try({
    model_output <- cj(data_subset, chosen_rec ~ Compensation,
                       id = ~ caseid,
                       weights = ~ weight,
                       estimate = "mm",
                       h0 = 0.5,
                       by = ~Energy)
    
    # Add group identifier
    model_output$threatccpol <- threatccpol_level
    
    # Bind results
    results_df <- bind_rows(results_df, model_output)
  }, silent = TRUE)
}


results_df2 <- results_df[!(results_df$threatccpol %in% c("Somewhat likely")),]


# Get unique values of threatccpol
threatccpol_vals <- unique(results_df2$threatccpol)

# Indexes where y-axis labels should be shown: 1st
show_indices <- c(1, 2)

# Create plot list
plot_list <- lapply(seq_along(threatccpol_vals), function(i) {
  o <- threatccpol_vals[i]
  show_y <- i %in% show_indices  # Only show y labels for specified positions
  
  ggplot(subset(results_df2, threatccpol == o), aes(x = estimate, y = level)) +
    geom_point(size = 3, color = "red") +
    geom_errorbarh(aes(xmin = lower, xmax = upper),
                   linewidth = 0.75, height = 0, na.rm = TRUE,
                   color = "red")  +
    # geom_text(aes(label = sprintf("%.2f", estimate)),
    #           vjust = -1.5, size = 3, color = "black") +  
    facet_wrap(~BY, nrow = 1) +
    scale_x_continuous(limits = c(0.3, 0.7), breaks = seq(0.3, 0.7, by = 0.1)) +
    geom_vline(xintercept = 0.5, linetype = "dashed", color = "gray50") +
    labs(
      subtitle = o,
      x = "Marginal Means",
      y = NULL
    ) +
    theme_bw(base_size = 12) +
    theme(
      plot.title = element_text(size = 10,face = "bold"),
      strip.text = element_text(size = 10),
      axis.text.x = element_text(size = 10),
      axis.title.x = element_text(size = 10),
      panel.grid = element_blank(),
      axis.text.y = if (show_y) element_text(size = 10) else element_blank(),
      axis.ticks.y = if (show_y) element_line() else element_blank()
    )
})

# Combine plots vertically
final_plot <- wrap_plots(plotlist = plot_list, ncol = 1) + 
  plot_annotation(
    title = "Interaction between energy and compensation by threat of CC policies"
  )
print(final_plot)

#Save plot
ggsave("figures/Figure_8.png",
       final_plot,
       width =8, height = 2.5*length(threatccpol_vals))

#------------------------------------------------------------------------------#
#### Environmental vulnerability ####
#------------------------------------------------------------------------------#

#Compute marginal means
affarea1 <- cj(
  long,
  chosen_rec ~ Energy*Compensation + Financing, id = ~caseid, by = ~affarea,
  weights = ~ weight, estimate = "mm"
)

#Define plot
affarea1plot <- plot(affarea1, group = "affarea", 
                     header_fmt = "%s",
                     vline = 0.5,
                     size = 3) +
  labs(title = "Support for policy package",
       subtitle = "by area of residence affected by pollution/extreme weather") +
  guides(colour=guide_legend(title="Affected by pollution/extreme weather"))  +
  aes(shape = feature) + # map group to `shape` aesthetic
  scale_shape_manual(values=c(15, 16, 17)) +
  ggplot2::theme(plot.title = element_text(size=20),
                 plot.subtitle = element_text(size=16),
                 axis.text=element_text(size=16),
                 axis.title=element_text(size=16),
                 text = element_text(size = 18), 
                 legend.position = "bottom",
                 legend.title = element_text(size=14),
                 legend.text = element_text(size=14),
                 axis.text.y = element_text(face=myFaces)
  ) +
  scale_shape(guide = 'none')


affarea1plot

#Save plot
ggsave(filename = 'figures/Figure_9.png',
       plot = affarea1plot, width = 12, height = 10)

#------------------------------------------------------------------------------#
#### Left-Right self-placement ####
#------------------------------------------------------------------------------#

#Compute marginal means
lr1 <- cj(
  subset(long, !is.na(long$lr6)),
  chosen_rec ~ Energy*Compensation + Financing, id = ~caseid, by = ~ lr6,
  weights = ~ weight, estimate = "mm"
)

#Define the plot
lr1plot <- plot(lr1[c(1:11, 23:33, 45:55),], group = "lr6",
                header_fmt = "%s",
                vline = 0.5, size = 3) +
  labs(title = "Support for policy package by Left-Right") +
  guides(colour=guide_legend(title="Ideology")) +
  aes(shape = feature) + # map group to `shape` aesthetic
  scale_shape_manual(values=c(15, 16, 17)) +
  ggplot2::theme(plot.title = element_text(size=20),
                 plot.subtitle = element_text(size=16),
                 axis.text=element_text(size=16),
                 axis.title=element_text(size=16),
                 text = element_text(size = 18), 
                 legend.position = "bottom",
                 legend.title = element_blank(),
                 legend.text = element_text(size=14),
                 axis.text.y = element_text(face=myFaces)
  ) +
  scale_shape(guide = 'none')


lr1plot

#Save the plot
ggsave(filename = 'figures/Figure_10.png',
       plot = lr1plot, width = 10, height = 10)

#------------------------------------------------------------------------------#
# Interaction between energy and compensation by left-right levels
#------------------------------------------------------------------------------#

# Initialize empty data frame
results_df <- data.frame()

#Loop over levels of lr6
for (lr6_level in unique(long$lr6)) {
  # Subset data
  data_subset <- subset(long, lr6 == lr6_level)
  
  # Skip if not enough data or no variation in variables
  if (nrow(data_subset) < 1 || length(unique(data_subset$Compensation)) < 2) next
  
  # Run model with error handling
  model_output <- try(
    cj(data_subset, chosen_rec ~ Compensation,
       id = ~ caseid,
       weights = ~ weight,
       estimate = "mm",
       h0 = 0.5,
       by = ~Energy),
    silent = TRUE
  )
  
  # Skip iteration if cj() failed
  if (inherits(model_output, "try-error")) next
  
  # Add lr6 label
  model_output$lr6 <- lr6_level
  
  # Bind results
  results_df <- bind_rows(results_df, model_output)
}

#Variable as factor for plotting
results_df$lr6 <- factor(results_df$lr6,
                         levels = c("Left", "Centre-left", "Centre", "Centre-right", "Right", "Not located"))

# Get unique values of lr6
lr6_vals <- levels(results_df$lr6)


# Indexes where y-axis labels should be shown: 1st, 3rd, 5th
show_indices <- c(1, 3, 5)

# Create plot list
plot_list <- lapply(seq_along(lr6_vals), function(i) {
  o <- lr6_vals[i]
  show_y <- i %in% show_indices  # Only show y labels for specified positions
  
  ggplot(subset(results_df, lr6 == o), aes(x = estimate, y = level)) +
    geom_point(size = 3, color = "red") +
    geom_errorbarh(aes(xmin = lower, xmax = upper),
                   linewidth = 0.75, height = 0, na.rm = TRUE,
                   color = "red") +
    # geom_text(aes(label = sprintf("%.2f", estimate)),
    #           vjust = -1.5, size = 3, color = "black") +  
    facet_wrap(~BY, nrow = 1) +
    scale_x_continuous(limits = c(0.25, 0.75), breaks = seq(0.25, 0.75, by = 0.25)) +
    geom_vline(xintercept = 0.5, linetype = "dashed", color = "gray50") +
    labs(
      title = o,
      x = "Marginal Means",
      y = NULL
    ) +
    theme_bw(base_size = 12) +
    theme(
      plot.title = element_text(size = 10,face = "bold"),
      strip.text = element_text(size = 10),
      axis.text.x = element_text(size = 10),
      axis.title.x = element_text(size = 10),
      panel.grid = element_blank(),
      axis.text.y = if (show_y) element_text(size = 10) else element_blank(),
      axis.ticks.y = if (show_y) element_line() else element_blank()
    )
})

# Combine plots vertically
final_plot <- wrap_plots(plotlist = plot_list, ncol = 2) + 
  plot_annotation(
    title = "Interaction between energy and compensation by Left-Right self-placement" 
  )
print(final_plot)

#Save plot
ggsave("figures/Figure_11.png",
       final_plot,
       width = 15, height = 2*length(lr6_vals))


#------------------------------------------------------------------------------#
# Predicted support for policy packages
#------------------------------------------------------------------------------#

#Create interaction between energy and compensation variables
long$Interaction1 <- interaction(long$Energy, long$Compensation)

#Compute marginal means for all combinations of factors
interfit3 <- cj(long, chosen_rec ~ Interaction1 ,
                id = ~ caseid,
                weights = ~ weight,
                estimate = "mm",
                h0 = 0.5,
                by = ~Financing)


#KEEP POLICY PARADIGMS SPECIFIED IN PAPER
packagedat <- interfit3 %>%
  filter(
    (BY == "Cuts in social spending" & level == "Rely on fossil fuels.Support affected companies") | #GROWTH FIRST
      (BY == "Cuts in social spending" & level == "Invest in renewables.Support affected companies") | #GREEN GROWTH
      (BY == "Wealth tax" & (level == "Invest in renewables.Support affected workers")) | #JUST TRANSITION
      (BY == "Wealth tax" & level == "Downscale production.Cash to all HHs") #POST GROWTH
  )

#Assign name to policy packages
packagedat$name <- NA


packagedat$name[packagedat$BY == "Cuts in social spending" & 
                  packagedat$level == "Rely on fossil fuels.Support affected companies"] <- "GROWTH FIRST\nFossil fuels\nSupport companies\nSocial cuts"


packagedat$name[packagedat$BY == "Cuts in social spending" & 
                  packagedat$level == "Invest in renewables.Support affected companies"] <- "GREEN GROWTH\nInvest in renewables\nSupport companies\nSocial cuts"


packagedat$name[packagedat$BY == "Wealth tax" & 
                  packagedat$level == "Invest in renewables.Support affected workers"] <- "JUST TRANSITION\nInvest in renewables\nSupport workers\nWealth tax"

packagedat$name[packagedat$BY == "Wealth tax" & 
                  packagedat$level == "Downscale production.Cash to all HHs"] <- "POST GROWTH\nDownscale production\nCash all HHs\nWealth tax"


#Transform "name" to factor variable
packagedat$name <- factor(packagedat$name,
                          levels = c("GROWTH FIRST\nFossil fuels\nSupport companies\nSocial cuts",
                                     "GREEN GROWTH\nInvest in renewables\nSupport companies\nSocial cuts",
                                     "JUST TRANSITION\nInvest in renewables\nSupport workers\nWealth tax",
                                     "POST GROWTH\nDownscale production\nCash all HHs\nWealth tax"))

table(packagedat$name)

#Define plot in Black and White
packageplot <- ggplot(packagedat, aes(x = name, y = estimate)) +
  geom_bar(stat = "identity", position = position_dodge(), colour = "black", fill = NA) +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2, position = position_dodge(0.9)) +
  labs(x = "Policy paradigm", y = "Marginal mean estimate") +
  theme_classic() +
  theme(legend.position = "none",
        text = element_text(size = 16), 
        axis.text = element_text(size = 12)) +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black")

packageplot

#Save the plot
ggsave(filename = 'figures/Figure_12.png',
       plot = packageplot, width = 10, height = 6)

#------------------------------------------------------------------------------#